
#include <navierstokes.hpp>


#include <feel/feelmesh/meshmover.hpp>
//#include <feel/feelalg/solvernonlinear.hpp>
//#include <feel/feelvf/exprbase.hpp>
//#include <feel/feelvf/expr.hpp>
#include <feel/feeldiscr/operatorinterpolation.hpp>

//#include <feel/feelalg/backendtrilinos.hpp>
//#include <Teuchos_ParameterList.hpp>

//#include <boost/mpl/count_if.hpp>
//#include <cldata.hpp>
//#include <appli_core.hpp>
//surtout pas #include <ns_cl.cpp>
#include <ns_cl.cpp>

#include <research/life-workspace/Fluid-Structure/functionSup.cpp>
#include <research/life-workspace/Fluid-Structure/methodsnum.hpp>

namespace Feel {
    namespace FSI {

        //---------------------------------------------------------------------------------------------------------//
        //---------------------------------------------------------------------------------------------------------//

        template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
        void
        NavierStokes<Dim, Order, Entity>::solveFluid(std::string __PDEtype)
        {

            if (M_LookComment) std::cout <<"[NavierStokes] : exec start \n";

            // vector solution of last step time
            vector_ptrtype Uvec( M_backend->newVector( M_Xh ) );
            *Uvec = *U_fluid;

            auto clDef = CL_DEF_NS(M_time);

            MethodsNum< self_ptrtype > meth(this->shared_from_this());

            //std::string PDEtype = "Stokes";
            if (__PDEtype == "Stokes")
                meth.linearSolver(clDef,Uvec);
            else if (__PDEtype == "Navier-Stokes")
                meth.AlgoNewton(clDef,Uvec);

            //meth.AlgoNewton2(clDef,Uvec);

            *U_fluid = *Uvec;

            if (M_LookComment) std::cout <<"[NavierStokes] : exec finsh \n";

        }

        //---------------------------------------------------------------------------------------------------------//
        template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
        void
        NavierStokes<Dim, Order, Entity>::updateLinearPDE(const vector_ptrtype& X,sparse_matrix_ptrtype& A , vector_ptrtype& F)
        {
            updateStokes(X,A,F);
        }

        //---------------------------------------------------------------------------------------------------------//

        template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
        void
        NavierStokes<Dim, Order, Entity>::updateALEmap()
        {

            if (M_LookComment) std::cout <<"[NavierStokes] : update ALE map start \n";

            using namespace Feel::vf;
            //structure_space_ptrtype Sh_top = M_depl_top->functionSpace();
            //structure_space_ptrtype Sh_bottom = M_depl_bottom->functionSpace();

            //Attention diferent celon que l'on donne le deplacement ou
            //la position de l'interface
            element_structure_type bc_disp_top( Sh_top, "bc_top" );
            //bc_disp_top = *M_depl_top;
            //bc_disp_top -= *M_bc_reference_top;
            bc_disp_top = *M_depl_top;
            bc_disp_top.updateGlobalValues();
            element_structure_type bc_disp_bottom( Sh_bottom, "bc_bottom" );
            bc_disp_bottom = *M_depl_bottom;
            //bc_disp_bottom -= *M_bc_reference_bottom;//vaut 0
            bc_disp_bottom.updateGlobalValues();
            std::vector<element_structure_type> polyBoundarySet;
            polyBoundarySet.push_back(bc_disp_bottom);
            polyBoundarySet.push_back(bc_disp_top);

            //define set of flags
            std::map< std::string, std::vector<flag_type> > flagSet;
            flagSet["fixed_bc"].push_back(mesh_fluid->markerName("Sortie"));//1
            flagSet["fixed_bc"].push_back(mesh_fluid->markerName("Entree"));//3
            flagSet["moving_bc"].push_back(mesh_fluid->markerName("ParoiB"));//2
            flagSet["moving_bc"].push_back(mesh_fluid->markerName("ParoiH"));//4

            // generate alemap with boundary on reference and actual mesh
            M_aleFactory->generateHighOrderMap( flagSet["moving_bc"], M_referencePolyBoundarySet, polyBoundarySet );

            *M_aux = M_aleFactory->getDisplacement();
            //essaye de comprendre pourquoi? !!!!!!!!!!!!!!!!!!!
            bool remesh = false;
            if ( remesh )
                {
                    mesh_fluid.reset(new mesh_type());
                    mesh_fluid = this->createFluidMesh( meshSize,"domainFluid" );
                }
            else
                {
                    *M_aux -= M_bdf_ale_map->unknown(0);
                }

            M_aux->updateGlobalValues();

            //deplacement du maillage
            M_mesh_mover.apply(mesh_fluid, *M_aux );

            if ( remesh )
                {
                    M_Xh = space_fluid_type::New( mesh_fluid );
                    U_fluid.reset(new element_fluid_type(M_Xh, "ub"));
                }

            M_identity_ale.reset(new element_fluid_velocity_type(M_Xh->template functionSpace<0>(), "identity_ale" ) );
            *M_identity_ale = vf::project( M_Xh->template functionSpace<0>(), elements( M_Xh->mesh() ), P() );
            M_identity_ale->updateGlobalValues();

            //std::cout << "\n vitesse mesh = " << M_bdf_ale_identity->polyDerivCoefficient(0) << "\n";

            *M_mesh_velocity = *M_identity_ale;
            //mesh_velocity.scale( M_bdf_ale_identity->polyDerivCoefficient(0) );
            //mesh_velocity.add( -1.0, M_bdf_ale_identity->polyDeriv() );
            *M_mesh_velocity -= M_bdf_ale_identity->unknown(0);
            M_mesh_velocity->scale( M_bdf_ale_identity->polyDerivCoefficient(0) );
            M_mesh_velocity->updateGlobalValues();

            //---------------------------------------------------------------

            if (M_LookComment) std::cout <<"[NavierStokes] : update ALE map finish \n";

        }

        //---------------------------------------------------------------------------------------------------------//
        //---------------------------------------------------------------------------------------------------------//

        template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
        void
        NavierStokes<Dim, Order, Entity>::init()
        {
            if (M_LookComment) std::cout<<"[NavierStokes] : init start\n";

            this->changeRepository();

            //deplacement initial
            *M_depl_top = vf::project( Sh_top, elements(Sh_top->mesh()), vf::cst(0.0) );
            *M_depl_bottom = vf::project( Sh_bottom, elements(Sh_bottom->mesh()), vf::cst(0.0) );

            if (M_cpt_export>0) {
                // m.a.j de la map ALE
                this->updateALEmap();
            }

            // initialisation des bdfs
            *M_identity_ale = vf::project( M_Xh->template functionSpace<0>(), elements( M_Xh->mesh() ), vf::P() );
            M_identity_ale->updateGlobalValues();

            M_bdf_ale_identity->shiftRight( *M_identity_ale );
            M_bdf_ale_map->shiftRight( M_aleFactory->getDisplacement() );
            M_bdf_fluid->shiftRight( *U_fluid );

            M_bdf_ale_map->initialize( M_aleFactory->getDisplacement() );
            M_bdf_ale_identity->initialize( *M_identity_ale );
            M_bdf_fluid->initialize( *U_fluid );

            M_bdf_ale_map->start();
            M_bdf_ale_identity->start();
            M_bdf_fluid->start();

            if (M_LookComment) std::cout<<"[NavierStokes] : init finish\n";

        }


        //---------------------------------------------------------------------------------------------------------//
        //---------------------------------------------------------------------------------------------------------//


        template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
        void
        NavierStokes<Dim, Order, Entity>::updateDeplacement(element_structure_ptrtype __depl_top,
                                                            element_structure_ptrtype __depl_bottom )
        {
            if (M_LookComment) std::cout<<"[NavierStokes] : updateDeplacement start\n";

            this->changeRepository();

            M_depl_top = __depl_top;
            M_depl_bottom = __depl_bottom;
            M_depl_top->updateGlobalValues();
            M_depl_bottom->updateGlobalValues();

            if (M_isMoveDomain)
                {
                    this->updateALEmap();
                }

            if (M_LookComment) std::cout<<"[NavierStokes] : updateDeplacement finish\n";
        }



        template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
        void
        NavierStokes<Dim, Order, Entity>::updateTime(double __t)
        {
            M_time = __t;
        }

        template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
        void
        NavierStokes<Dim, Order, Entity>::updateBdf()
        {
            if (M_LookComment) std::cout <<"[NavierStokes] : updateBDF start\n";

            boost::timer btime;
            btime.restart();

            this->application()->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                                                   % this->application()->about().appName()
                                                   % convex_type::name()
                                                   % Order
                                                   % this->application()->vm()["fluide-hsize"].template as<double>()
                                                   );

            M_bdf_ale_identity->shiftRight( *M_identity_ale );
            //M_bdf_ale_map->shiftRight( /*M_aleFactory->getDisplacement()*/*M_aux );
            M_bdf_ale_map->shiftRight( M_aleFactory->getDisplacement() );
            M_bdf_fluid->shiftRight( *U_fluid );

            M_bdf_ale_identity->next();
            M_bdf_ale_map->next();
            M_bdf_fluid->next();

            double timeElaps = btime.elapsed();
            if (M_LookComment) std::cout << "[NavierStokes] : updateBDF finish in " << timeElaps << "s\n";
        }


        template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
        typename NavierStokes<Dim, Order, Entity>::element_fluid_ptrtype
        NavierStokes<Dim, Order, Entity>::getNormalStress()
        {

            using namespace Feel::vf;

            if (M_LookComment) std::cout << "[NavierStokes] : getNormalStress start\n";

            auto solFluid = getSolution();

            element_fluid_0_type u = solFluid->template element<0>();
            element_fluid_1_type p = solFluid->template element<1>();

            //Tenseur des deformations
            auto defv = 0.5*(gradv(u)+trans(gradv(u)) );
            //Matrice Identitee
            auto Id = oneX()*trans(oneX()) + oneY()*trans(oneY());
            // Tenseur des contraintes (trial)
            auto Sigmav = -idv(p)*Id + 2*M_mu*defv;

            //element_fluid_type Stress(M_Xh,"Stress");
            element_fluid_ptrtype Stress(new element_fluid_type (M_Xh,"Stress"));
            element_fluid_0_type uStress = Stress->template element<0>();
            modifVec222(uStress, Stress, Sigmav*N() , "ParoiB" );
            modifVec222(uStress, Stress, Sigmav*N() , "ParoiH" );

            M_Ufluid_ref->container() = Stress->container();

            if (M_LookComment) std::cout << "[NavierStokes] : getNormalStress finish\n";

            return M_Ufluid_ref;//->template element<0>();
        }


        template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
        void
        NavierStokes<Dim, Order, Entity>::runPtFixe()
        {

            using namespace Feel::vf;

            std::cout << "ptFixe start\n";

            //----------------------------------------------------------------------------------------------------------//

            if ( this->application()->vm().count( "help" ) )
                {
                    std::cout << this->application()->optionsDescription() << "\n";
                    return;
                }

            //----------------------------------------------------------------------------------------------------------//

#if 0
            Teuchos::ParameterList My_List;
            My_List.set("Max Iters", 1550);
            My_List.set("Tolerance", 1e-10);
            My_List.set("Solver", "GMRES");

            Teuchos::ParameterList& Prec_List = My_List.sublist("Preconditioner");
            Prec_List.set("Type", "ILU");
            Prec_List.set("Drop Tolerance", 1e-3);
            cout<< "\nOLA" <<My_List << endl;

            //( static_cast<BackendTrilinos>(*M_backend)).set_options(My_List);
            BackendTrilinos bbb;
            bbb.set_options(My_List);
#endif

            double erreur=1;
            double err_tol=1e-6;//this->vm()["tolerP"].template as<double>();
            double err1,err2,norm1,norm2;
            uint NbMaxItPtFixe=30;
            uint cptItPtFixe=0;
            uint cptTime=0;


            //U_fluid->container() = M_bdf_fluid->unknown(0).container();

            element_fluid_ptrtype Uoldb( new element_fluid_type(M_Xh,"uoldb"));
            erreur=1;
            cptItPtFixe=0;

            vector_ptrtype R( M_backend->newVector( M_Xh ) );
            vector_ptrtype T( M_backend->newVector( M_Xh ) );
            vector_ptrtype S( M_backend->newVector( M_Xh ) );


            vector_ptrtype F( M_backend->newVector( M_Xh ) );
            sparse_matrix_ptrtype A( M_backend->newMatrix(M_Xh,M_Xh) );

            element_fluid_velocity_type uYYY(M_Xh->template functionSpace<0>(),"xxx");

            while (erreur>=err_tol && cptItPtFixe < NbMaxItPtFixe)
                {
                    *Uoldb = *U_fluid;
                    //std::cout<<"\nLe Uold Norm L2 : " << Uoldb->l2Norm() << "\n";
                    //this->updatePtFixe(*U_fluid,*Uoldb,*M_mesh_velocity,*M_aux);

                    this->updatePtFixe(*Uoldb,A,F);
                    //this->updateStokes(*Uoldb,A,F);

                    std::cout << "PtFixe : Solve start\n";
                    vector_ptrtype X( M_backend->newVector( M_Xh ) );
                    M_backend->solve( _matrix=A, _solution=X, _rhs=F/*,_maxit=23456*/ );
                    //bbb.solve( A,A,X,F/*,_maxit=23456*/ );
                    *U_fluid = *X;
                    std::cout << "PtFixe : Solve finish\n";


                    *T = *U_fluid;
                    //this->updateResidual( T, R );
                    //modifVec222(uYYY,R, vec(cst(0.),cst(0.)),"Entree");
                    //modifVec222(uYYY,R, vec(cst(0.),cst(0.)),"ParoiH");
                    //modifVec222(uYYY,R, vec(cst(0.),cst(0.)),"ParoiB");
                    //std::cout<<"\nLe residu du ptfixe Norm L2 : " << R->l2Norm() << "\n";

                    erreur = /*R->l2Norm();*/this->diffL2(*U_fluid,*Uoldb);
                    std::cout << "Erreur : " << erreur <<"\n";
                    ++cptItPtFixe;

                }
            std::cout << "ptFixe finish\n";
        }

    } // end namespace NavierStokes
} // end namespace Feel


// instantiation
template class Feel::FSI::NavierStokes<2,2,Feel::Simplex>;
//template class Feel::NavierStokes<2,4,Feel::Simplex>;


